Electromagnetic modes of a disordered photonic crystal 
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We present a systematic numerical approach to compute the eigenmodes and the related eigen- 
frequencies of a disordered photonic crystal, characterized by small fluctuations of the otherwise 
periodic dielectric profile. The field eigenmodes are expanded on the basis of Bloch modes of the 
corresponding periodic structure, and the resulting eigenvalue problem is diagonalized on a trun- 
cated basis including a finite number of Bloch bands. The Bloch-mode expansion is very effective for 
modeling modes of very extended disordered structures in a given frequency range, as only spectrally 
close bands must be included in the basis set. The convergence can be easily verified by repeating 
the diagonalization for an increased band set. As illustrations, we apply the method to the Wl 
line-defect waveguide and to the L3 coupled- cavity waveguide, both based on a photonic crystal 
slab with a triangular lattice of circular holes. We compute and characterize the eigenfrequencies, 
spatial field profiles and radiation loss rates of the localized modes induced by disorder. For the 
Wl waveguide, we demonstrate that radiation losses, at the bottom of the spatially-even guided 
band, are determined by a small hybridization with Bloch modes of the spatially-odd band, induced 
by disorder in spite of their frequency separation. The Bloch-mode expansion method has a very 
broad range of applications: it can be also used to accurately compute the modes of structures with 
systematic modulations of the periodic dielectric constant, as in several designs of high-Q cavities. 

PACS numbers: 42.70.Qs, 42.25.Dd, 72.15.Rn, 78.20.Bh 



I. INTRODUCTION 

Thanks to the periodic spatial dependence of their di- 
electric constant, photonic crystal (PHC) structures 1 3 
make it possible to engineer the confinement and prop- 
agation of the electromagnetic field, and hold great 
promise for applications to integrated optics. 

Similarly to atomic crystals, PHCs always present 
some amount of disorder as a result of the fabrication 
process. 4 When considering PHC slabs - a periodic array 
of cylindric holes in a dielectric slab waveguide - disorder 
arises from the irregular shape of the holes and from the 
roughness of the dielectric interfaces. Disorder is mostly 
considered as a disturbance, in the context of applica- 
tions. It prevents ballistic light propagation, by inducing 
losses due to scattering, 5 16 and possibly strong local- 
ization of light, 5 ' 16 24 thus constituting a major physi- 
cal limit to the use of PHCs as interconnects in pho- 
tonic circuits, and in all applications relying on slow light 
propagation. 25 27 Disorder-induced scattering is also re- 
sponsible for extrinsic radiation losses, 6,28 ' 29 as it mixes 
Bloch momenta of truly guided and leaky modes of the 
nominally periodic structure. 

More generally, it has been suggested that PHCs might 
be an election system for studying the fundamental prop- 
erties of light scattering in disordered media. It is rather 
singular that the seminal work by Sayeed John, 2 consid- 
ered as one of the first works to introduce the idea of 
a PHC, actually considers the photonic bands only as a 
tool to achieve Anderson localization of light. 30 The idea 
put forward in the work by John is that, in the vicinity 
of a point of high symmetry, the dispersion of a photonic 
band can be approximated by a parabola and Maxwell 
equations can consequently be mapped onto an effective 



Schrodinger equation in a continuous effective medium, 
including a random potential term derived from the ac- 
tual disorder profile of the system. 

Modeling the effects of disorder in PHCs is a current 
research topic and a very challenging computational task. 
In a regular structure, Bloch's theorem allows restricting 
the solution of Maxwell equations to a single elemen- 
tary cell. Taking advantage of the crystal periodicity, 
several computational approaches have been developed, 
based on mode expansion 28,31 ' 32 or scattering matrix 33 36 
formalisms. Disorder however breaks the translational 
symmetry of the crystal, and Bloch's theorem can no 
longer be applied. Then, a full solution of Maxwell equa- 
tions over an extended region of space, possibly includ- 
ing thousands of elementary cells would be needed both 
for studying fundamental properties, such as Anderson 
localization, and for modeling device operation. First 
principle calculations based on a finite-element sampling 
of real space are thus of very limited use, as the computa- 
tional cost of simulating such extended samples would be 
prohibitive. Theoretical approaches to disordered PHCs 
have been proposed, based on various levels of approx- 
imation. Some properties, such as diffusive light prop- 
agation and the corresponding radiation losses, can be 
predicted already within a perturbation theory of light 
scattering. 6 ' 16 ' 29 ' 37 More generally, a variety of coupled- 
mode approaches have been proposed for studying dis- 
ordered PHCs. 8 ' 12 ' 14 ' 38 A coupled-mode theory in gen- 
eral accounts for light scattering between few selected 
modes. The reduced number of modes makes it pos- 
sible to model the scattering processes in a fully non- 
perturbative fashion, but restricts the predictive power 
of the method. As an example, Anderson localization 
of light arises from multiple coherent scattering between 
many different propagating modes, 22 ' 30 ' 39 and typically 
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requires a solution of Maxwell equations beyond few- 
mode coupling, to be modeled correctly. It should be 
noted however that the idea of mode coupling can be 
generalized to scattering theories including a full basis of 
propagating modes. In this case, the theory is formally 
equivalent to a full solution of Maxwell equations, as is 
e.g. the case for the quasi- normal Bloch-mode scattering 
matrix approach recently developed by Lecamp et al. 9 ' 40 

A general aspect of most theoretical approaches 33 36 ' 40 
is the fact that they compute directly the linear response 
of the periodic dielectric medium. For several reasons 
however, it would be beneficial to compute the actual 
Maxwell eigenmodes. Apart from the direct visualization 
of the Anderson-localized spatial profiles, knowledge of 
the actual eigenmodes and of their detailed Bloch-mode 
components can give insight into the mechanisms of slow- 
light propagation, 25 ' 27 ' 41 44 radiation 6 ' 28 ' 29 and disorder- 
induced 5 16 losses. When modeling systems with coupled 
photonic and electronic degrees of freedom, such as quan- 
tum dots embedded in cavities or guides, 19 ' 45 53 or PHC 
polaritons, 31 ' 54 it is most natural to start with the eigen- 
states of both coupled subsystems, especially within a 
fully quantum mechanical treatment where second quan- 
tization of electromagnetic modes is needed. If the linear 
response function is known, then eigenmodes can be ob- 
tained by finding the poles of its analytical continuation 
on the complex frequency plane. 35 Conversely, if Maxwell 
eigenmodes are known, the linear response function can 
always be expressed as a function of these eigenmodes, by 
using the resolvent representation in Fredholm theory. 55 
Both indirect ways are however computationally far from 
being optimal. A method for directly computing Maxwell 
eigenmodes of an extended PHC with a small perturba- 
tion of the periodic dielectric profile would thus be an 
indispensable complement to linear response theories. 

Here, we develop the Bloch-mode expansion formalism, 
for directly computing Maxwell eigenmodes of a photonic 
structure characterized by a small perturbation of the 
otherwise periodic dielectric profile. The basic idea be- 
hind this formalism is that eigenmodes of the perturbed 
structure can be expanded on the basis of Bloch modes of 
the regular, unperturbed PHC. This expansion extends, 
in general, to all bands and all Bloch momenta of the 
regular structure. If the perturbation is small, as in the 
case of disorder in state-of-the-art PHC slabs, then the 
expansion can be truncated on a relatively small subset of 
Bloch bands - spectrally close to the frequency region of 
interest - and the Maxwell equation can be diagonalized 
restricting to this subspace, at a considerably reduced 
computational cost. The convergence of the method can 
be easily assessed by repeating the calculation on a pro- 
gressively augmented Bloch-band subset. After laying 
down the general formalism, as examples of application 
we study the Wl line-defect waveguide 56 ' 57 and the L3 
coupled-cavity waveguide, 58 both based on a PHC slab 
structure with a triangular lattice of circular holes. For 
computing the Bloch modes of the regular PHC, we use 
the guided-mode expansion method recently proposed by 



Andreani and Gerace, 28 because it offers a good bal- 
ance between predictivity and ease of implementation, 
although some limitations in its accuracy to predict out- 
of-plane losses have been recently put forth. 9 This is not 
however the only possible choice, and the Bloch-mode 
expansion could equally rely on a first-principle numer- 
ical solution of Maxwell equations for computing Bloch 
modes of the regular structure. We show how a fully 
converged calculation of modes in the guided region of 
both waveguides can be achieved at a minimal compu- 
tational cost, for structures as long as 1024 elementary 
PHC cells. We study the localization of the modes and 
their radiation loss rates, and discuss the influence of dis- 
order and losses on slow-light propagation, in the light of 
the existing literature. For the Wl waveguide, we show 
that radiation losses at the bottom edge of the spatially- 
even guided band are determined by hybridization with 
the spatially-odd guided band, in spite of their mutual 
frequency separation. This result explains the physical 
limitations to the Q-factor of Wl-based structures, such 
as "gentle confinement" cavities, 59 61 and shows the way 
to a significant optimization of these systems. 

In Section II, we present the theoretical formalism. 
We also discuss the method for computing radiation loss 
rates, when Bloch modes are computed via guided-mode 
expansion. In Sections III and IV, we discuss the results 
of the application of the present method to a Wl line- 
defect waveguide and to a L3 coupled-cavity waveguide, 
respectively. Section V presents the main conclusions 
and an outlook of future applications. 



II. THEORY 

We consider a PHC made of a periodic array of cylin- 
drical air holes etched in a dielectric slab of thickness d. 
The regular PHC - namely in the absence of disorder - 
is described by the dielectric constant e(r). In presence 
of disorder instead, the dielectric constant is defined as 

C '(r) = 6(r) + <fe(r) (1) 

As the holes are present only in the slab, we assume that 
both e(r) and e'(r) take constant values in the upper and 
lower semi-infinite claddings. In the following, we will 
always consider an "air bridge" configuration, for which 
e = e' = 1 and Se = in the claddings. Here, we are 
neglecting the roughness that might characterize the slab 
interfaces, as we are primarily interested in the disorder 
affecting the cylindrical holes. 

The starting point of our analysis are the Bloch modes 
Ekn(r) of the regular PHC structure under study. They 
are the eigenmodes of the Maxwell equation 

V x V x E kn (r) - 4^(r)E kn (r) = , (2) 

where uj^ n are the corresponding eigenfrequencies. These 
modes might have been obtained by any computational 
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method that takes advantage of Bloch's theorem. For ex- 
ample, they could be computed by exact finite-element 
solution of Maxwell equations, or by an approximate 
scheme like e.g. guided- mode expansion. 28 In the ex- 
amples of application that follow, we will consider a 
finite system size with periodic boundary conditions, 
thus implying a discrete set of Bloch-momentum values 
k. The formalism presented in this Section can how- 
ever be straightforwardly generalized to a continuous k- 
spectrum. 

In presence of disorder, the actual eigenmodes of the 
system are denoted by E^(r), and are indexed by a global 
index /?, as the Bloch momentum k is no longer a con- 
served quantity. These modes are solution of the Maxwell 
equation 



UJ 



V x V x E^r) - ^(rjE^r) = . (3) 



They can be formally expanded on the basis of the Bloch 
modes of the regular structure: 



E^r) = ^^(k,n)E k „(r). 



(4) 



kn 



By replacing (4) into (3), and using (1) and (2), we obtain 
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-£e(r) - J-5e(r) 



E kn (r) =0. 



We take the scalar product of Eq. (5) by E£, n/ (r) and 
integrate on r (we assume the modes Ek n (r) to be nor- 
malized over the total volume of the system). By taking 
advantage of the orthogonality relation 62 

J dre(r)E k , n ,(r) • E k „(r) = <W<W , (6) 

we finally obtain 



\a'n' _ 



UJ- 



kn 



-UJ}. 



^(k',n') = 0. 



(7) 

Eq. (7) is a generalized eigenvalue problem, formally 
equivalent to the full solution of Maxwell equation (3). 
The matrix elements Vk n ^> n > contain all information 
about the disorder profile. They are defined as 



= J dr<5e(r)E k , n ,(r) • E k „(r) . 



(8) 



It should be pointed out that, when assuming periodic 
boundary conditions for the Bloch modes, the Bloch- 
mode expansion (4) implies periodic boundary conditions 
also for the modes E^(r). Thus, boundary effects at the 
side terminations of the PHC (e.g. interface with a ridge- 
waveguide), that might give rise to additional spectral 
features in some structures, 10 should be computed addi- 
tionally, e.g. through mode-coupling or transfer-matrix 
approaches that go beyond the scope of the present work. 



It is interesting to recall the principle underlying the 
present formalism had partly been laid down already in 
the seminal work by S. John. 2 There, it is pointed out 
that, close to a band edge, the Maxwell problem can be 
approximated by an effective Schrodinger equation whose 
kinetic term is essentially determined by the curvature of 
the single Bloch band under study. The potential term is 
nonlocal and is given by the matrix element of the ran- 
dom potential between actual electric field modes. This 
exactly corresponds to our Eqs. (7) and (8) if restricted 
to the edge of a single Bloch band. According to the basic 
theory of electronic states in crystals, restricting to the 
edge of a single Bloch band coincides with the effective- 
mass approximation. The present formalism generalizes 
John's original idea to an expansion over several bands. 
Its advantage lies in the possibility to truncate the vector 
space defining the eigenvalue problem (7) to a subset of 
the bands of the regular PHC, in order to ensure optimal 
convergence of the properties under study. In particu- 
lar, when studying the effect of disorder on the band no, 
it should be possible to neglect all bands n whose fre- 
quencies are much further away from that of band no, 
than the disorder contribution (V). By inspection of Eq. 
(7), it is easily shown that bands can be neglected if the 
condition 



^k n 



•^knl > ^knl^k no,kn| 



(9) 



is fulfilled. The reliability of the method can be assessed 
by repeating the diagonalization of (7) for an increas- 
ing number of bands and checking for convergence of the 
eigenvalues and eigenvectors of interest. This perturba- 
tion criterion holds however only for n ^ no. Within 
the same band - in particular for low- dimensional sys- 
tems - even the slightest amount of disorder can induce 
a strong mixing of Bloch modes in correspondence to 
band extrema, possibly giving rise to Anderson local- 
ization, as originally suggested by John. 2 The typical 
situation in which condition (9) is particularly favor- 
able is that of defect states in the bandgap of a two- 
dimensional PHC, for which most other bands are nat- 
urally separated by the frequency bandgap. Among the 
systems falling in this class are all one-dimensional sys- 
tems based on a two-dimensional PHC, such as line- 
defect waveguides, 56 ' 57 coupled-cavity waveguides, 63 ' 64 
etc. These systems are of particular interest because of 
the potentially high quality factor and of the possibil- 
ity of slow-light propagation. 25 ' 27 Slow-light in particu- 
lar can be exploited to enhance the coupling of light to 
the electronic states of embedded quantum dots, for ap- 
plications as single photon emitters. 19 ' 45 53 In this con- 
text in particular, knowing the actual eigenmodes of 
the electromagnetic field is crucial for modeling their 
coupling to electronic states. 55 Otherwise, the method 
should also be useful to model the states of a disordered 
two-dimensional PHC at the edges of the bandgap. 17 ' 20 
Finally, it is also possible to apply the present method 
to special designs involving a small systematic perturba- 
tion Se(r) of a periodic PHC, as in the case of long defect 
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cavities 60 or cavities obtained by a slow grading of the 
lattice parameters. 59 61 In the next section we will apply 
the present method to a Wl line-defect waveguide in a 
PHC made of a triangular lattice of air holes, and to the 
related system of a L3 coupled-cavity waveguide. In the 
case of the Wl waveguide, we will show that accounting 
for more than one Bloch band is essential in determining 
the properties of the Maxwell eigenmodes. 

One of the main consequences of disorder in PHCs 
are the extrinsic radiation losses. In a regular struc- 
ture, modes Ek n (r) lying outside the light cone in the 
(c^k, k)-space are truly guided within the slab and ide- 
ally lossless. Inside the light cone, on the other hand, 
the electromagnetic modes are characterized by a contin- 
uous frequency spectrum, with a spectral function that 
has narrow resonances in correspondence to the band dis- 
persion cjkn- In this case, formally the index n only la- 
bels the resonance frequencies, while eigenmodes should 
be labelled by a continuous frequency variable. The con- 
tinuous frequency spectrum is related to the lossy na- 
ture of these modes, and implies a finite rate of radiation 
out of the plane into the claddings. When disorder is 
present, Eq. (4) mixes discrete and continuous parts of 
the spectrum, resulting in a continuous spectrum over 
the whole (cjk? k)-space. Now all modes (4) have compo- 
nents in the continuous part of the spectrum, inside the 
light cone, and are therefore lossy. For these lossy modes, 
Eqs. (4)- (8) would not strictly hold and should be mod- 
ified to account for the continuous frequency spectrum. 
Then, a formally exact but computationally demanding 
approach to evaluate the loss rates, would require solv- 
ing the Maxwell eigenvalue problem (2) numerically in 
three dimensions, and investigating the properties of its 
continuous spectrum. Here we prefer to adopt a compu- 
tationally less demanding approach, based on the guided 
mode expansion recently proposed by Andreani et al. 28 
to compute the modes of the regular structure Ek n (r). 
These modes are expanded on a subset of the full Maxwell 
basis, formed by the guided modes of an effective homo- 
geneous dielectric slab 

E k „(r) = c»(k + G, a)E^ G » , (10) 

where a runs over the different guided modes. The leaky 
modes of the effective slab are thus neglected in this ex- 
pansion. This approximation turns out to be highly pre- 
dictive in typical situations, at a marginal computational 
cost. Radiation losses of Bloch modes can then be evalu- 
ated by a perturbative approach analogous to the Fermi 
golde rule in quantum mechanics, involving overlap inte- 
grals between modes (10) and leaky modes of the effective 
slab. The continuous spectrum of leaky Bloch modes is 
thus accounted for by an irreversible decay of the approx- 
imately guided modes (10) into the leaky modes of the 
homogeneous slab. The same approach can be applied in 
the present context to compute the radiation loss rates of 
the eigenmodes E^(r) of the disordered structure. Then, 
to the real eigenvalues oo 2 o/c 2 is associated an imaginary 



part T/3 related to the loss rate. The expression for this 
quantity is 

k,G A=TE,TM j=l,3 V / 

where g = k + G, pj is the density of leaky modes for 
a given frequency and in-plane momentum, and is 
the overlap integral between mode E^(r) and the leaky 
mode of the effective slab with momentum g. The index 
j runs over the two cladding layers, while the index A 
runs over the polarization basis of the radiative modes. 
From here, the actual loss rate can be computed as 

lp * c 2 ry (2 W/3 ) . (12) 

The explicit expressions for the quantities entering Eq. 
(11) are presented in the Appendix. One limit of validity 
of this perturbative approach is of course that 7^ <C cup. 

In presence of disorder, an additional condition must 
be fulfilled. While in a regular structure light would 
propagate ballistically in each Bloch mode Ek n , disor- 
der is generally responsible for damping of the trans- 
mission along the plane, induced by multiple scattering 
between different Bloch momenta and in particular co- 
herent backscattering. 9 ' 12 This damping is typically de- 
scribed as a disorder-induced loss mechanism. The same 
multiple scattering can induce Anderson localization of 
modes E^(r), and the average attenuation length of the 
transmission at a given frequency is a measure of the av- 
erage localization length of modes at that frequency. In 
this picture, a competition between radiation losses and 
multiple scattering arises. The destructive interference 
in multiple scattering, at the origin of forward attenua- 
tion and Anderson localization, requires that radiation 
losses are negligible over the total multiple scattering 
path. As we plan to compute modes E^(r) based on 
a guided mode expansion of Bloch modes Ek n (r) - thus 
without accounting for radiation losses - then an addi- 
tional condition of validity of the present approach is 
that the attenuation length corresponding to 7/3 be much 
longer than the localization length of modes E^(r). If 
this is not the case, then Eq. (4) would be inconsistent 
with the radiation loss rate 7/3 computed perturbatively 
afterwards, as the Bloch waves involved in the expan- 
sion would be damped before producing the interference 
pattern that gives rise to the localized mode. This condi- 
tion holds particularly well for modes that in the regular 
structure would lie outside the light cone and thus be ra- 
diation lossless. A typical example are the guided modes 
of a Wl line-defect waveguide, as discussed in the next 
Section. 

Before proceeding to applications, we point out that 
the present method is fully non-perturbative. If the 
Bloch-mode basis, on which the expansion is made, is 
computed exactly, then the converged expansion (4) coin- 
cides with the full solution of the Maxwell problem for the 
extended disordered structure. If, as we do in the exam- 
ples below, the Bloch basis is computed via guided-mode 
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expansion, 28 then modes (4) are still close to the full so- 
lution of Maxwell equations, with the only assumption 
that the admixture of out-of-plane propagating modes is 
negligible and can be taken into account perturbatively 
at a later stage through Eq. (11). 



III. APPLICATION TO A Wl LINE-DEFECT 
WAVEGUIDE 

As an example of application of the present method, 
we consider the problem of disorder in Wl line-defect 
waveguides. We choose this system, as it was widely 
investigated both theoretically and experimentally, thus 
representing a valid platform on which our method can 
be tested. The effect of disorder on one-dimensional sys- 
tems is particularly dramatic, as all states are expected 
to be exponentially localized. 39 In the context of PHC 
waveguides, the role of disorder was highlighted already 
in early experimental studies. Two disorder-related ef- 
fects are of primary importance and were extensively 
discussed in the literature. First, the restrictions on 
slow-light propagation, due to coherent multiple scat- 
tering that induces an attenuation of the transmitted 
intensity. 5 16,25 ' 27 ' 41 44 Second, the appearance of extrin- 
sic radiation losses. 6 ' 28 ' 29 In what follows, we characterize 
the eigenmodes of a Wl waveguide in presence of disor- 
der and discuss their properties in the light of the existing 
studies. We show that the present method allows to ac- 
tually compute the localized eigenmodes, to study their 
localization length and their radiation loss rates, and to 
characterize the disorder-induced transmission losses. In 
addition, the method provides new insight into the physi- 
cal mechanisms governing radiation losses. In particular, 
we show that extrinsic radiation losses in Wl waveguides 
are mostly determined by hybridization of the spatially- 
even and spatially-odd guided Bloch modes - a fact that 
has great impact on the design and optimization of high- 
Q optical cavities based on "gentle confinement" of Wl- 
guided modes. 59 61 

We consider a Wl line-defect waveguide, consisting in 
one missing row of holes in a triangular lattice of cylindri- 
cal air holes of radius r, with lattice parameter a and slab 
thickness d. The Wl waveguide is sketched in Fig. 1(a). 
The modes Ek n (r) of the regular structure are computed 
by means of guided- mode expansion. 28 The waveguide is 
oriented along the x-axis, while the z-axis is perpendicu- 
lar to the slab. We adopt a rectangular supercell of width 
a along x and height ba\^3/2 along y, where b is an in- 
teger, as sketched in Fig. 1(a). Modes that are strongly 
confined into the guide core, are scarcely affected by the 
superperiodicity along y introduced by the supercell. 56 
We can thus safely restrict to modes having k y = 0, cor- 
responding to an effective one-dimensional system whose 
bands are characterized by Bloch momentum k x along 
the x-axis. For simplicity, in what follows we denote 
k x simply as k. For the guided- mode expansion, we re- 
strict to the lowest guided mode of the effective slab, 
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FIG. 1: (a) Sketch of the Wl line-defect waveguide. The 
shaded region denotes the supercell used for computing the 
bands of the regular structure, (b) Sketch of the disorder 
model used in this work. The dashed and full circles indicate 
holes of the regular and disordered structure respectively, (c) 
Band structure of the Wl waveguide with d = 0.5a, e = 12, 
r = 0.3, and b = 10. Band 1 (red) is the spatially-even guided 
band. Bands labeled 1 to 4 are those that were included in 
the Bloch-mode expansion of the disordered structure. 



having TE polarization. We thus neglect the TE-TM 
mixing induced by the hole lattice - an approximation 
that holds very well for small d/a, when the frequency 
of the next guided mode falls outside the band gap of 
the 2-D lattice. 28 Calculations can be easily generalized 
to an expansion on several guided modes if needed, by 
a straightforward extension of the formalism. As usual, 
a scaling relation with respect to the lattice parameter 
a holds and we can express all lengths, momenta and 
frequencies in dimensionless units. For the calculations 
that follow, we set d/a = 0.5, r/a — 0.3, b — 10, and 
we consider an air-bridge structure with dielectric con- 
stant ei 5 3 = 1 in the cladding layers, and 62 = 12 for the 
slab material, corresponding to an effective-slab dielec- 
tric constant 62 = 8.77 for the calculation of the guided 
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modes (see Appendix). In the calculations we use Ho's 
method 56 for the Fourier transform of the dielectric pro- 
file, and we truncate the expansion at G ma x = 3 (in 
units of 27r/a), corresponding to a basis of 229 guided 
modes. Fig. 1(c) shows the band structure in the re- 
gion of the guided mode. Bands labeled 1 to 4 will be 
included in the truncated band set used for the disorder 
calculations. 80 Band 1 is the main guided band, often 
referred to as index-guided or spatially-even band (with 
respect to the &k z mirror plane), while band 2 is the gap- 
guided, spatially-odd band. 

We model disorder as a gaussian fluctuation of both 
hole radii r m and positions (x m ,y m ), where m runs over 
all holes of the extended structure, as sketched in Fig. 
1(b). In particular, r m = r + oY m , x m = x$ + 6x m , and 

Vm = Vm +Sy m , where r and (x$,y$) are respectively 
the radius and the positions of the holes in the regu- 
lar structure. The fluctuations 5r m , 5x m , and Sy m are 
Gauss-distributed with standard deviations <r r , <j x , and 
(j y . This "minimal" disorder model - in particular the 
radii fluctuations - is expected to catch the main effects 
of hole shape fluctuations on the spectral properties 4 ' 9 ' 10 , 
at a low numerical cost in computing the disorder matrix 
elements Vkn,k'n', as shown in the Appendix. More re- 
fined disorder models, possibly accounting for the rough- 
ness of the hole contour 6 ' 12 ' 65 ' 66 , could be easily imple- 
mented. As discussed above, the converged expansion (4) 
corresponds to the exact solution of the Maxwell problem 
and would thus automatically include local field effects 
in proximity of the hole sidewalls, whose importance has 
recently been debated. 31 ' 65 ' 66 In what follows, we restrict 
our analysis to the situation in which a x = a y — a r = a. 
Although we didn't carry out a systematic study, our sim- 
ulations seem to confirm that radii fluctuations are the 
dominant effect, while position fluctuations give a less 
significant contribution. We consider values of a ranging 
from 0.001a to 0.008a, corresponding to state-of-the-art 
structures. 4 ' 10 

For the numerical simulations, we consider Wl struc- 
tures of length L = 1024a, namely including 9216 holes 
for the super cell height b = 10 considered here. This cor- 
responds to assuming a uniform grid of 1024 ^-points in 
the interval [— 7r/a,7r/a] for the eigenvalue problem (7). 
A single disorder realization of a system of this length 
requires a few hours on a desktop computer. Computing 
Maxwell eigenmodes of such an extended structure with 
other first-principle methods, like e.g. a finite-element 
sampling of real space, would constitute a practically 
untractable problem. The advantage of the present ap- 
proach is therefore striking. 

In order to study the spectral properties of the system, 
a natural choice would be to introduce the quantity 



\Hf}(k,ky =0)| S 

U - LOj3 + ijp 



(13) 



/3-th eigenmode, and 7/3 is the rate of radiation losses 
computed according to (11). Eq. (13) is a scalar Green's 
function of the electromagnetic field, corresponding to 
the trace of the Green's tensor of Maxwell equations. 67 It 
is evaluated at equal initial and final Bloch momenta (for- 
ward propagation) , thus providing information about the 
transmission amplitude and, ultimately, about the spec- 
tral density of the field. 39 It is expressed through the 
resolvent representation 67 based on the magnetic field 
modes H#(r). These modes constitute a complete or- 
thonormal basis with the standard measure for the scalar 
product, as opposed to the dielectric measure e(r) needed 
to define a scalar product for electric field modes. 55 ' 62 
They are related to the corresponding electric field modes 
asH^(r) = — i(c/^)VxE^(r). In Eq. (13) however, the 
choice k y = corresponds to integrating the field along 
y, thus missing all spectral features related to spatially- 
odd modes with respect to the &k z mirror symmetry. In 
order to be able to visualize these spectral contributions, 
we prefer to adopt a different choice and introduce the 
quantity 



Gy(k,U)) 



|Hg(fc,y)| 2 

^ UJ - UJ/3 + 27/3 



(14) 



where Hp(k,k y = 0) is the Fourier transform, taken at 
k y = (and z = 0) of the magnetic field H^(r) of the 



where Hp(k,y) is obtained by Fourier transforming the 
field H#(r) with respect to x only. Eq. (14) represents 
the Green's function for the propagation of the field along 
x, at a fixed position y. In this way, by evaluating (14) at 
2/^0, i.e. displaced from the symmetry plane &kz, spec- 
tral signatures of both even and odd modes will appear. 
For all plots, we choose y = 0.067a. 

Fig. 2(a) is a log-scale intensity plot of the spectral 
density A(k,uj) = -Im(G y (fc,cj)) for a = 0.002a, ob- 
tained by diagonalizing Eq. (7) on the subspace spanned 
by bands 1 and 2. The dispersion curves of the two bands 
for the regular structure are superimposed, as well as 
the free-space light dispersion. This latter defines the 
boundary of the light cone: below it, modes of the di- 
electric slab are fully confined by total internal reflection 
and ideally lossless. The spectral width along band 1 
shows an abrupt increase above the light cone, due to 
the onset of the intrinsic radiation loss mechanism. A 
similar abrupt change is visible along band 2 at about 
ka/27r — 0.4. On the left of that point on the plot, the 
frequency of band 2 becomes degenerate with points of 
the same band lying inside the light cone, thus maximiz- 
ing the mode mixing due to disorder and the subsequent 
extrinsic radiation losses, as already suggested by Ku- 
ramochi et al. 7 As a consequence, the further increase of 
the spectral width across the light cone boundary is less 
abrupt for this band, but still noticeable. Figs. 2(b)-(e) 
are linear intensity plots of the same quantity in the vicin- 
ity of the bottom of the guided band 1, for a = 0.001a, 
0.002a, 0.004a, and 0.008a, respectively. A feature com- 
mon to the four cases is the presence of sharp spectral 
signatures below the bottom of the band of the regular 
structure. These features are extended along and ex- 
tend to increasingly lower frequencies as a is increased. 
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FIG. 2: Spectral function (see text) of the disordered Wl 
waveguide, (a) a = 0.002a, plotted on a log 10 color scale. 
Bands 1 and 2 of the regular structure are also plotted (full 
lines). The dotted line is the boundary of the light cone, (b)- 
(e) Detail (on a linear color scale) of the spectral function 
at the lower edge of band 1, for a = 0.001a, a = 0.002a, 
a = 0.004a, a = 0.008a, respectively. 



They originate from spatially localized modes induced 
by disorder, and contribute to form the well known Lif- 
shitz tail in the density of states below the band of the 
regular crystal. 68 70 These features have been experimen- 
tally characterized by Le Thomas et al. 42 and the corre- 
sponding sharp features in the transmission spectra have 
been experimentally observed and theoretically studied 
by modeling the optical response. 10,12 

The density of states can be computed by averaging 
over several statistical realizations of the disorder config- 
uration. Fig. 3(a) shows the density of states obtained 
by averaging over 500 realizations, for a = 0.002a, com- 
pared to the density of states of the regular structure 
(panel (b) shows the spectral density for reference). The 
Van Hove singularity of the regular crystal is smeared 
over an interval of the order of uocr /(2irc) and the Lifshitz 
tail is clearly visible. 

The main innovation brought by the present approach 
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FIG. 3: (a) Density of states in the vicinity of the edge of band 
1. Blue: a = 0.002a, averaged over 500 disorder realizations. 
Red: the corresponding quantity for the regular structure, 
showing the characteristic Van Hove singularity at the band 
edge, (b) The spectral function and the band dispersion are 
plotted as a reference. 



is the possibility to compute the field profile of the ac- 
tual Maxwell eigenmodes. In connection to the study of 
disorder, this makes it possible to study the localization 
properties and brings insight to the question whether An- 
derson localization of light can occur in PHCs. 2 ' 10 ' 18 24 
We plot the electric field profile of three selected modes 
of a disorder realization with a = 0.002a in Fig. 4(a)- (c) 
respectively. The three selected frequencies are indicated 
by dotted lines in Fig. 2(c). The main plot in each 
panel represents the value of the field intensity at the 
center of each elementary cell of the guide (in this way 
the main envelope can be plotted over the full length of 
the simulated structure), while the insets show the full 
electric field profile for selected representative regions of 
the structure. Fig. 4(a) shows the ground eigenmode at 
ua/27rc = 0.27267, which lies well below the band edge 
and is localized over a few elementary cells of the struc- 
ture. Fig. 4(b) shows a mode at ua/27rc = 0.27283, lying 
very close to the band edge and displaying a less local- 
ized character with few main lobes. Most of the modes 
computed in this spectral region have a similar character, 
with several lobes and an overall envelope with exponen- 
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FIG. 4: (a)-(c) Electric field spatial profile of three selected 
eigenmodes computed for a = 0.002a, at uja/2-KC — 0.27267, 
0.27283, and 0.27518, respectively (indicated by horizontal 
lines in Fig. 2(c)). The main panels display the field inten- 
sity computed at the center of each elementary cell, while 
insets show the full field profile for selected portions of the 
Wl waveguide. 



tially decaying tails, as expected for Anderson localized 
states. Modes with several lobes have been studied in 
various kinds of PHCs and are considered as precursors of 
the extended necklace states. 10 ' 23 ' 24 ' 42 ' 71 ' 72 Finally, Fig. 
4(c) shows a mode at ujaj2i\c = 0.27518, well within the 
energy range of the guided band, which is delocalized 
over the whole length of the simulated sample while still 
showing irregular amplitude fluctuations induced by dis- 
order. This is expected as in the spectral region of this 
mode, the localization length is much longer than the 
simulation length L considered here. 

The radiation loss rate 7/3 for each eigenmode j3 has 
been computed using Eqs. (11) and (12). They are 
plotted, for one disorder realization, in Fig. 5(a), for 
a = 0.001a, 0.002a, 0.004a, and 0.008a. Panel (b) in this 
figure shows the band structure and the light line for 
reference. Let us focus first on frequencies below 0.287, 
for which the guided band of the regular structure lies 
outside the light cone. For this frequency range, radia- 
tion losses are extrinsic, as already discussed. For each 
given value of the disorder amplitude a, the rates 7/3 
behave quite regularly, except for the lowest frequencies 
in the Lifshitz tail of the spectrum. There, the values 
of 7/3 are significantly scattered over about one decade, 
as one would expect for modes that are strongly local- 
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FIG. 5: (a) Radiation loss rates of disordered eigenmodes, 
computed for the four values of a considered in this work. 
Modes for which the intrinsic radiation loss mechanism is 
dominant, are those showing a loss rate independent of a. 
(b) Band dispersion (full) and light cone (dotted) are plotted 
for reference. 



ized by the random fluctuations of the underlying crys- 
tal structure. A similar trend, with strong variation of 
the loss rates close to the band edge, was obtained in 
Ref. 29, although the perturbative approach to disor- 
der adopted in that work cannot predict the fluctuations 
of 7/3. These large fluctuations are typical of states in 
the localized part of the spectrum, namely where the 
localization length is much smaller than the simulation 
length L. Similar fluctuations are observed in the trans- 
mission spectrum 9 ' 10,12 and are more generally expected 
for all properties of localized eigenmodes, as in the local- 
ization lengths discussed below. Above the fluctuating 
region, the rates 7/3 are approximately constant over a 
wide frequency range, until they dramatically increase 
when approaching the boundary of the light cone. The 
rates 7/3, for the values of the disorder amplitude con- 
sidered here, are in general very small, as expected for 
an extrinsic process. When the disorder amplitude a is 
doubled, the extrinsic rates correspondingly increase by 
a factor of four. This dependence has already been dis- 
cussed in the literature. 11 It can be explained by consid- 
ering that the overlap integral entering Eq. (11) depends 
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FIG. 6: (a)-(c) Eigenfrequencies of four selected eigenmodes 
computed as a function of the number of bands included in 
the Bloch-mode expansion, for a = 0.002a. (d) Computed 
radiation loss rates of the same modes, as a function of the 
number of included bands, (e) Coefficients of the Bloch-mode 
expansion \Up(k, n)\ computed for the fundamental mode f3 = 
1 and n = 1,2,3,4. 



linearly on the disorder matrix element (V) to leading 
order. This can be inferred from Eq. (A9), if we consider 
a perturbation expansion of mode where the zeroth 
order term is a guided Bloch mode, that does not con- 
tribute to the overlap integral (A9) due to orthogonality. 
For the frequency range above the light line, the eigen- 
modes of the disordered structure are characterized by 
much higher loss rates that are now mostly of intrinsic 
nature and no longer depend on a. 81 In the frequency 
range between 0.294 and 0.306, rates of modes related 
to the spatially-odd band 2 are also present. Here losses 
show again a partially extrinsic character, especially in 
the lowest frequency part where a quadratic dependence 
on (7 appears. For the discussion below, it is also impor- 
tant to keep in mind that rates in this frequency range are 
significantly higher than those computed at the bottom 
of band 1. 

The origin of the radiation loss rates can be inferred 
from the nature of the computed eigenmodes. These, 
according to Eq. (4) are linear combinations of Bloch 
modes at different k and different bands. This hybridiza- 
tion involves all modes of the regular structure, regard- 
less of their frequency spacing Au - the mixing fraction 
being quantified by (V)(uj/Auj) to leading order. Fig. 
6(a)-(c) shows the convergence of the eigenfrequencies of 
four selected eigenmodes, at frequencies corresponding to 



band 1, as a function of the number of bands included 
in the diagonalization. The convergence is excellent al- 
ready when only band 1 is included in the Bloch-mode 
expansion. Convergence of the loss rates, shown in Fig. 
6(d), requires instead that the spatially-odd band 2 is 
also included in the calculation. In particular, for the 
three lowest frequencies lying in the region of extrinsic 
radiation losses, the loss rates increase by as much as 
two decades when band 2 is included. The reason is that 
the main source of extrinsic losses, close to the edge of 
band 1, is the hybridization with Bloch modes of the 
spatially-odd band 2. This mixing would be forbidden 
by symmetry in the regular structure and is now allowed 
as disorder breaks the a^z invar iance. The fourth selected 
frequency instead corresponds to an eigenmode for which 
the intrinsic radiation loss mechanism is dominant, and 
the hybridization with band 2 doesn't thus produce any 
significant effect. Fig. 6(e) shows the components of the 
lowest computed eigenmode /3 = 1, \Up(k,n)\. Apart 
from the dominant contribution of Bloch modes close to 
k = ±0.5 of band 1, the main contribution comes from 
band 2 that shows about 3 x 10 -3 admixture at the edge 
of the Brillouin zone. This explains the dominant con- 
tribution of band 2 to the extrinsic loss rates, seen in 
Fig. 6(d). Band 2 losses are several orders of magnitude 
higher than the extrinsic losses computed when account- 
ing for band 1 only. Therefore, even the slightest admix- 
ture of the two bands - here about 0.3% - can completely 
determine the final loss rate at the band edge. The role 
of the mixing between spatially-odd and even bands has 
been already pointed out by Kuramochi et al. 7 in the 
framework of perturbation theory, thus restricted to the 
frequency region where the two modes are degenerate. 
Our result generalizes the one by Kuramochi et al., by 
showing that the mixing determines all extrinsic radia- 
tion loss rates of the main guided band of a Wl waveg- 
uide. The importance of this finding goes well beyond 
the study of disorder effects. It suggests that radiation 
losses in all Wl-based systems can be reduced by study- 
ing designs that minimize the mixing between bands 1 
and 2, e.g. by displacing band 2 at higher frequencies. 
It can indeed be argued that the same mechanism limits 
the Q-factor of cavities based on the gentle confinement 
of the Wl guided mode. 59 61 Our conclusion thus shows a 
new direction to take for the optimization of such cavity 
designs. 

To conclude the analysis of the Wl waveguide, we 
study the localization lengths of the simulated eigen- 
modes. The localization length in a disordered system 
can be defined in several different ways. 39 One natural 
definition, when studying eigenstates, is the inverse par- 
ticipation number (IPN). For a wavefunction (j)(x) in 1-D, 
the IPN is defined as I = 1/ J \(f)(x)\ 4 dx , and has the di- 
mension of a length if | </>(#) | 2 * s normalized to 1. The 
IPN measures the cumulative length over which the field 
intensity is significantly large. For very singular shapes 
of the disorder profile, the IPN can differ considerably 
from the exponential decay length of the tails of <p{x) - 
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FIG. 7: (a) Circles: IPN values of the lowest eigenmodes com- 
puted for one disorder realization with a = 0.002a. Blue line: 
IPN averaged over 500 disorder realizations. Red line: Av- 
erage attenuation length of the transmission along the guide. 
Inset: detail in the vicinity of the edge of band 1 (denoted 
by a vertical line), (b) The spectral function and the band 
dispersion are plotted as a reference. 



which is another measure of the localization length. 39 For 
an eigenmode H#(r) of our system, we thus define 



dx\Hp(x,y = 0)\ 4 



(15) 



where, in order to extract the mode length, we take the 
integral along the guide core at y = 0. A second char- 
acteristic length related to disorder, typically studied in 
the context of PHCs, 8 10 ' 12 ' 22,23 is the exponential at- 
tenuation length of the transmission of a monochromatic 
beam. For this, we introduce the quantity 



T(x, x' , uj) 



U*Jx,y = O)-U (x',y = O) 



LO - Up + H/3 



(16) 



which is proportional to the transmission of a monochro- 
matic source from x to x' . The frequency-dependent 
exponential attenuation length is defined as the decay 
constant L u in T(x,x',uj) oc exp( — \x r — x\/L UJ ). In the 
numerical calculation of both L u and Ip , we checked the 



convergenge of the computed values as a function of the 
system size L, up to the largest value of L considered, in 
order to detect the onset of finite-size effects. Fig. 7(a) 
shows the quantities Ip and L^, as computed for the Wl 
guide under study. Lines denote configuration averages 
over 500 disorder realizations, while the circles are IPN 
values for the eigenmodes of a given realization. The 
IPN is more affected by the finite size of the simulation 
and saturates at about Ip = 250a, while the saturation 
of the attenuation length is not visible on the scale of 
the plot. This is expected, as the IPN is an extensive 
property of the eigenmode profile, contrarily to the at- 
tenuation length which can be extracted from the slope 
at a single point. When finite-size effects are negligible, 
the two quantities roughly coincide. We point out that 
the attenuation length includes the decay due to radi- 
ation losses, while the IPN is an exclusive property of 
the eigenmodes of Eq. (7). In the frequency range of 
Fig. 7, the computed losses correspond to a negligible 
contribution to the attenuation length of the transmis- 
sion, which is thus uniquely determined by the coherent 
multiple scattering in the propagation along the guide. 
As for other properties of localized eigenmodes, the IPN 
displays large fluctuations at each given frequency. This 
suggests that, independently of the average localization 
length at a given frequency, rare necklace modes can exist 
with a profile that extends over the whole sample length 

£ 10,23,24,42,71,72 

Localization in one-dimensional photonic 
structures has been reported in several 
experiments. 10 ' 18 ' 19 ' 21 ' 23 ' 24 ' 42 The present result sheds 
light on the question concerning the nature of light 
localization in these structures. We have already 
excluded the extrinsic radiation losses as a cause of 
attenuation in a Wl waveguide close to the bottom 
of band 1. It is clear that the origin of localization 
in our model calculations is to be found in coherent 
multiple scattering processes, corresponding to the 
Bloch-mode superposition in expansion (4). This is 
of course expected, as in one-dimensional disordered 
systems localization is always present, 30 ' 39 and the Wl 
waveguide shows to be very closely one-dimensional. 
Indeed, the expansion on a single 1-D Bloch band was 
sufficient to obtain well converged eigenmodes, while the 
very small admixture with band 2 was only necessary to 
account correctly for radiation losses. A relevant ques- 
tion, raised by Vlasov et al., 22 is whether localization is 
of the Anderson kind or rather due to a random local 
energy shift of the band edge induced by large disorder 
fluctuations. This second mechanism can be pictured 
as the random alternance in space of pass-band and 
stop-band regions, and is thus enforced by the Bragg 
diffraction of the underlying periodic structure, rather 
than by scattering in an effective medium as in the 
Anderson localization case. Our approach gives direct 
access to the spectral function. It shows that, for the 
disorder amplitudes considered here and corresponding 
to state-of-the-art samples, 4 ' 10 the band-edge shift is 
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much smaller than the width of the stop band (see Figs. 
1 and 2). According to Vlasov et al. then, at least for 
the disorder model considered here localization should 
be of the Anderson kind. Comparing the behavior of the 
IPN and the attenuation length below the band edge, as 
shown in the inset of Fig. 7(a), supports this conclusion. 
For decreasing frequency, the attenuation length ceases 
to follow the average IPN and decreases down to its 
smallest possible value, essentially given by the unit 
cell length a. For frequencies below this crossover, 
the attenuation length is no longer determined by the 
localization length of the eigenmodes, but rather by 
Bragg diffraction on the PHC, and a similar attenuation 
length would arise within the stop band also in the case 
of a regular crystal. The IPN, on the other hand, seems 
to approach a constant value of a few unit cells, that 
corresponds to the smallest IPN taken by Anderson 
localized modes for the given value of <r, as already 
observed when directly visualizing the mode in Fig. 
4(a). 



IV. APPLICATION TO A L3 
COUPLED-CAVITY WAVEGUIDE 

As an additional test of our method, we consider a 
chain of L3 coupled cavities, as sketched in Fig. 8(a). 
The system differs from the Wl waveguide by the pres- 
ence of one hole every four along the waveguide. This 
particular kind of coupled- cavity waveguide has recently 
been experimentally characterized, 58 in order to deter- 
mine the origin of out-of-plane radiation losses and the 
limiting factors to slow light propagation. More gen- 
erally, coupled-cavity waveguides have been the object 
of various experimental 58,73,74 and theoretical 13,63,64,75,76 
studies, because of their potential as effective slow-light 
channels, thanks to the very flat band expected for an 
impurity chain. 

Here, we start our analysis from the Bloch modes of 
the regular structure, that we compute as for the Wl 
waveguide using guided- mode expansion 28 , by assuming 
the supercell sketched in Fig. 8(a). The disorder model 
and all other parameters are the same as for the Wl 
waveguide studied above. In particular, the same length 
L = 1024a, corresponding to 256 L3 cavities, is assumed. 
With these parameters, the band of interest 58 goes from 
ua/27rc = 0.2766 to ua/27rc = 0.2779, as plotted in Fig. 
8(b). Calculations using our method fully converge al- 
ready when including only this band in the Bloch-mode 
expansion. This holds also for the radiation loss rates, 
contrarily to the case of the Wl waveguide, because of the 
dominant role of the intrinsic radiation loss mechanism 
as discussed below. The spectral function for a = 0.002a 
is plotted in Fig. 8(b), while Fig. 8(c) displays the 
radiation loss rates computed for the field eigenmodes, 
for the disorder amplitude varying from a = 0.001a to 
a = 0.008a. The Brillouin zone, four times smaller than 
for a Wl waveguide, now lies completely above the light 
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FIG. 8: (a) Sketch of the L3 coupled cavity waveguide. The 
shaded region denotes the supercell used for computing the 
bands of the regular structure, (b) Spectral function com- 
puted for a = 0.002a, showed as a linear scale color plot. 
The dispersion of the guided band of the regular structure is 
plotted as a white line, (c) Radiation loss rates of disordered 
eigenmodes, computed for the four values of a considered in 
this work. 



cone, and radiation losses are mostly determined by the 
intrinsic mechanism predicted for the regular system. 13,75 
In particular, the rates vary from a value lower than that 
of the single L3-cavity, 28,77,78 at the band bottom, to a 
higher one at the band top, as predicted in the absence 
of disorder from the 3-D solution of Maxwell equations 13 
or from a tight-binding approach. 75 . The dependence 
on the disorder amplitude a clearly shows that, up to 
a = 0.004a, the extrinsic radiation loss mechanism takes 
over the intrinsic one only in the vicinity of the band 
bottom. For larger <j instead, the extrinsic mechanism 
becomes relevant over the whole band, and large fluctu- 
ations of the rates appear at each given frequency. 

Fig. 9 shows the electric field intensity computed for 
three selected eigenmodes, denoted by horizontal lines in 
Fig. 8(b). The ground eigenmode in Fig. 9(a) is localized 
over the length of a few L3-cavities. Modes at higher fre- 
quencies displayed in panels (b) and (c) are progressively 
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FIG. 9: (a)-(c) Electric field spatial profile of three selected 
eigenmodes computed for a = 0.002a. The three correspond- 
ing values of uja/2itc are denoted by horizontal lines in Fig. 
8(b). The main panels display the field intensity computed 
at the center of each elementary cell, while insets show the 
full field profile for selected portions of the L3 coupled-cavity 
waveguide. 



more delocalized, again as expected for a one-dimensional 
disordered system. Notice that each mode has the shape 
of a slowly varying envelope function modulating the field 
profile of the ground mode of a L3 cavity, 77 as in tight- 
binding calculations. 75 

As for the Wl waveguide studied in the previous Sec- 
tion, the modes displayed in Fig. 9 are computed with- 
out accounting for losses [which are computed at a later 
stage through Eq. (11)]. Contrarily to the Wl waveg- 
uide however, here the attenuation length corresponding 
to the computed radiation loss rate is actually shorter 
than the localization length of the modes, when consid- 
ering the least localized modes at the band center. Hence, 
for these modes the result of our Bloch-mode expansion 
is not fully consistent with the initial assumption of loss- 
less Bloch modes. To better understand this issue, we 
show in Fig. 10(a) the localization lengths computed 
from the IPN and the attenuation length computed from 
the decay of the transmission (16). Fig. 10(b) shows the 
spectral function for reference. The attenuation length 
is evaluated both by accounting for the computed loss 
rates 7/9 (black curve) and by assuming a vanishing loss 
rate for all modes (red curve). This latter coincides with 
the average IPN over a wide range of frequencies close 
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FIG. 10: (a) Circles: IPN values of the eigenmodes computed 
for one disorder realization with a — 0.002a. Blue line: IPN 
averaged over 500 disorder realizations. Red line: Average 
attenuation length of the transmission along the guide. Inset: 
detail in the vicinity of the lower band edge (denoted by a ver- 
tical line) . (b) The spectral function and the band dispersion 
are plotted as a reference. 



to the band extrema, while at the band center the IPN 
starts being affected by the finite length of the simu- 
lated structure, similarly to the case of the Wl waveg- 
uide. The attenuation length including losses is instead 
dramatically suppressed, except for a narrow region at 
the band bottom (see inset). We argue that, in the re- 
gion where the two curves differ significantly, the local- 
ized eigenmodes computed here are actually ill-defined. 
Radiation loss rates of Bloch modes are so high that they 
should be accounted for already in the multiple scatter- 
ing process from which the eigenmodes of the disordered 
system arise. This could be achieved by replacing the 
guided-mode expansion with a full expansion on guided 
plus radiative modes of the slab, or simply with an expan- 
sion on 3-D plane waves. The inset of Fig. 10(a) shows a 
detail of the region below the band edge, displaying the 
same trend already observed for a Wl waveguide in Fig. 
7(a), although now the minimal IPN value is larger than 
in the Wl case, as expected from a longer unit cell. 

We conclude by making some considerations on the 
influence of disorder on slow-light propagation, both in 
the Wl and in the coupled-cavity waveguide. Slow-light 
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FIG. 11: Double logarithmic plot of the attenuation length 
of the transmission as a function of the nominal group index. 
The horizontal line indicates the limit set by the finite simu- 
lation length, while the diagonal line is proportional to n~ 2 
and is plotted as a guide to the eye. Blue: Wl waveguide. 
Red: L3 coupled-cavity waveguide neglecting radiation losses. 
Green: L3 coupled-cavity waveguide accounting for radiation 
losses. 



has been one of the main driving forces of the develop- 
ment of one-dimensional PHC structures. 25 ' 27 ' 44 ' 79 Ide- 
ally, close enough to the bottom of the guided band, the 
group velocity of a narrow light pulse, defined as v g = 
duJk/dk, can become arbitrarily small. Disorder hinders 
this ideal behavior by producing a strong attenuation of 
the transmitted signal, as discussed above. Past exper- 
imental studies 5 ' 7 ' 12 ' 42 have focused on the dependence 
of the attenuation length on the group index, defined 
as rig = c/vg. All recent theoretical studies 6 ' 9 ' 10 ' 12 ' 14 ' 16 
show that the dominant functional dependence of the at- 
tenuation length close to the band bottom is on n~ 2 . 
This behaviour is induced by disorder through the mech- 
anism of coherent backscattering, determined by second 
and higher-order coherent scattering processes, which is 
a precursor of localization. 39 The n~ 2 -dependence was 
confirmed in several experiments. 5 ' 7 ' 11 It was also pointed 
out 9 ' 10 ' 16 that, very close to the band bottom, the onset 
of localization makes multiple scattering beyond second 
order important and the n~ 2 -dependence consequently 
breaks down for high n g , with the attenuation length 
that levels off to its minimal value given by the lattice 
parameter a. In Fig. 11, we plot the computed atten- 



uation length L u for a = 0.002a as a function of the 
nominal group index n gi both for the Wl and for the 
coupled cavity waveguides. By "nominal", it is meant 
the group index that would result at frequency uo from 
the dispersion of the ideal band. The plot thus does 
not include frequencies below the band bottom. The re- 
sult for the Wl waveguide (blue) reproduces the result 
obtained by Mazoyer et al. 9 ' 10 For the coupled cavities, 
we display two curves corresponding to the case includ- 
ing radiation losses (green), and to the unphysical but 
illustrative lossless situation (red) respectively. This lat- 
ter coincides with the Wl result, except for the largest 
values of n g considered, that arise close to the band bot- 
tom. The interpretation of this result is straightforward. 
A small group index corresponds to a frequency well 
above the band bottom, where eigenmodes have a lo- 
calization length much longer than the elementary cell. 
In this case, independently of the length of the elemen- 
tary cell, the envelope function of the eigenmodes will be 
localized according to the John mechanism, 2 with the 
localization length solely determined by the curvature 
of the Bloch-mode dispersion, thus by the group index. 
When approaching the band bottom instead, the local- 
ization length becomes smaller and eventually reaches 
its lower bound corresponding to a small number of el- 
ementary cells, as shown from the analysis of the IPN. 
In this limit, the localization length of the coupled-L3 
system is four times longer than that of the Wl waveg- 
uide, at corresponding n g . This limiting behavior can 
be inferred by comparing the upper branch of the red 
curve (corresponding to the band bottom) and the blue 
curve in Fig. 11. When accounting for intrinsic radiation 
losses in the L3-coupled system, the attenuation length 
is dramatically reduced except for the region close to the 
band bottom, where losses are the smallest and disor- 
der concurs significantly in determining the attenuation. 
Following Mazoyer et al., 9 ' 10 , a diagram of the kind dis- 
played in Fig. 11 sets the maximum length scale over 
which light propagation can still be ballistic, and thus 
the slowdown predicted for the ideal system can still oc- 
cur. Then, from this analysis we can conclude that, in 
spite of the increased loss rates, a L3-coupled waveguide 
can still bring an advantage over the Wl waveguide when 
close to the band bottom, by increasing four times this 
maximum length for a given (large) group index. 



V. CONCLUSIONS 

The main result of the present analysis is the demon- 
stration that Boch-mode expansion can be an extremely 
effective and reliable tool for studying the properties of 
the eigenmodes of Maxwell equations, in PHCs where 
disorder or small systematic variations of the otherwise 
periodic structure are present. 

The approach relies on a method to compute the Bloch 
modes of the regular PHC. Ideally, this first elemen- 
tary step should be accomplished by exact diagonaliza- 
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tion of the Maxwell problem within a single elementary 
cell, for example by using a finite element approach. For 
the scope of the present work however, we preferred to 
adopt the guided- mode expansion method, 28 because of 
its reliability at a minimal computational cost. In this 
case, however, out-of-plane propagating modes are not 
included in the expansion, and the radiation losses of the 
final Maxwell eigenmodes must be computed non self- 
consistently within a perturbative approximation at a 
later stage. 

The general advantage of Bloch-mode expansion lies 
in the optimal choice of the Bloch-mode basis, for situ- 
ations where the variations from the periodic structure 
are small. Then, an expansion restricted to a very lim- 
ited number of spectrally close Bloch bands, can pro- 
duce excellent convergence, thus practically reproducing 
the exact solution of Maxwell equations at a reasonable 
computational cost. 

As an example, we have applied the method to study 
the effect of disorder in a Wl line-defect waveguide and 
a L3 coupled-cavity waveguide. The method gives direct 
access to the shape and spectrum of the eigenmodes, from 
which the localization properties, and the mechanism un- 
derlying radiation losses have been characterized. 

One very important result of this work - apart from the 
assessment of the Bloch-mode expansion method - is the 
explanation of the origin of radiation losses of the spa- 
tially even guided modes of a Wl waveguide close to the 
band bottom. Radiation losses are mainly determined by 
the hybridization of these modes with those of the guided 
band of opposite spatial parity. A very small admixture 
of the two bands is sufficient to increase the radiation loss 
rate close to the band edge by two decades. Hence, the 
effect is dominant in spite of the frequency separation of 
the two bands that hybridize - a feature overlooked by 
previous analyses. This result was only possible thanks 
to the separate analysis of the contribution of different 
Bloch bands. It unveils the physical mechanism limiting 
the quality factor of all structures based on a Wl waveg- 
uide. In particular, all designs of high-Q cavities based on 
slow modulation of a Wl waveguide should be affected. 
We suggest that a significant increase in the quality factor 
of such structures can be obtained by studying a waveg- 
uide design that reduces the hybridization, e.g. by maxi- 
mizing the frequency separation between the two guided 
bands. 

The Bloch-mode expansion has a very broad spectrum 
of potential applications. These range from the realistic 
modeling of high-Q cavities, to the study of radiation- 
matter coupling and cavity quantum electrodynamics ef- 
fects for systems of several quantum dots coupled to the 
modes of PHC cavities or waveguides. All these appli- 
cation should take advantage of the possibility of com- 
puting the actual eigenmodes of Maxwell equations, thus 
making the Bloch-mode expansion an election method 
for the modeling and design of photonic systems. 
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Appendix A: Guided-mode expansion 

The modes Ek n (r) of the regular PHC, and their eigen- 
frequencies uo\^ n can be computed using the guided-mode 
expansion method introduced by Andreani and Gerace 28 . 
The modes are then expanded on the basis of the guided 
modes of an effective homogeneous dielectric slab, char- 
acterized by dielectric constants efj, with j = 1, 2, 3 indi- 
cating the lower cladding, the PHC layer and the upper 
cladding respectively. The quantities ej are the spatial 
averages of the actual dielectric profiles in each layer. 
Guided modes E^ G a (r) are labeled by their total in- 
plane momentum g = k+G and mode index a = 1, 2, . . .. 
Then, the PHC modes are expanded as 

Ekn(r) = £c„(k + G,a)E<fj GiQ (r). (Al) 

G,a 

In this work we give explicit expressions only for guided 
modes having TE polarization, as we assume quasi-TE 
PHC modes in the results that we present in Sections III 
and IV. More general expressions, for modes where TE 
and TM polarizations are mixed, can be obtained in a 
straightforward way by following the guidelines in Ref. 
28. 

For the calculation of the matrix elements Vk^k'n'? 
only the expression of Ek n (r) for r in the slab is needed, 
as we assumed Se(r) = in the claddings. By further 
assuming a structure symmetric with respect to a mirror 
reflection a xy , namely equal upper and lower claddings, 
the guided-mode expansion reads 

Ekn(r) = 2^c n 0)i-^e g — 7 ^2Re(A 2fI ) cos(q^z) , 
/i c 

(A2) 

where /i = (k + G, a) is the global index for the guided 
mode, p is the in-plane position vector and z the coordi- 
nate orthogonal to the plane, = ^J^oj^/c 2 — g 2 , and 

S is a normalization area that will disappear from the 
final result. The complex coefficient Ai^ is defined in 
Ref. 28. We have introduced the unit polarization vec- 
tor e g = z x g of the guided mode, where z is the unit 
vector along the z-axis. By replacing (A2) into (8), we 
obtain 

W'n' = ^<(/i)c n VMg' - g)^W , (A3) 
/x,/x' 
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where 

M m , = 2^^e s -e g ,[A^A2^l2-+Re(A 2ft A2^)Ii^4) 

J 2± = 2 ™[(fr*v)l] . (A5) 

(Ifl ± V 

In Eq. (A3) enters the Fourier transform of the dielectric 
perturbation, defined as 

Se(g) = | j dpj*"6e{p,z = 0) . (A6) 
J s 

Here, we make use of the assumption, introduced in Sec- 
tion III, that the disorder profile Se(r) is constant along z 
within the slab. Within the disorder model introduced in 
that Section, circular holes of varying radius and position 
are assumed, and the Fourier transform (A6) can be eval- 
uated analytically. We start from the Fourier transform 
of the disordered dielectric profile e'(p, 0) 

e'(g) = I /d P e^V(p,0) 

= e 2 6 s , + (1 - e 2 ) ge^-^MgrJ^r) 

m 

where r m and p m = (x m ,y m ) are respectively the radius 
and the in-plane position of the ra-th hole of the struc- 
ture, Ji is the Bessel function of the first kind of order 
1, and the sum runs over all holes present in the PHC. 
An analogous expression holds for the Fourier transform 
of the regular dielectric profile e(g), if r m and p m are 

replaced by r and the positions p$ = (x$ , yffl ) of the 
holes in the regular structure. Then, Eq. (A6) is simply 
expressed as 

Se(g) = e'(g) - e(g) . (A8) 

We point out that the present method can be at least 
formally extended to more elaborate disorder models, 



with holes of non-circular shape and possibly with a z- 
dependence of the dielectric profile (e.g. tapered holes 
or roughness along z). For this, it is enough to account 
for the corresponding Se(r) in the integral (8). For an 
arbitrary disorder shape however, this integral must in 
general be computed in a fully numerical fashion, which 
might make the method computationally very demand- 
ing. 

In order to compute the radiation loss rate of each 
mode, we generalize the expression introduced by An- 
dreani et al. 28 . In analogy to the Fermi golden rule in 
quantum mechanics, here the loss rate is determined by 
its overlap matrix element with leaky modes of the ef- 
fective slab, through Eq. (11). Within the guided- mode 
expansion, the matrix element Wj^ reads 

= J / <Mr)E£(r) • E$(r) , (A9) 

where E^-(r) are the leaky modes of the effective slab, 
described in detail in Ref. 28, A runs over the two po- 
larizations TE and TM, and j = 1, 3 runs over the lower 
and upper cladding layer. By using the two expansions 
(4) and (Al), we can express the matrix element as 

m d = EE w> n )<( k ' + G '> a )K't > ( Ai °) 

k',nG' 5 a 

where 

2 

K>t = ^ / dre(r)E$» • E$(r) . (All) 

Because of the translational symmetry of the lattice, the 
above integral is proportional to #kk' , thus lifting the sum 
over k 7 in (A10). Details for the evaluation of Eq. (All) 
are given in Ref. 28. 
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Here, as in Ref. 29, discontinuities appear, due to the sin- 
gularities in the 1-D density of states entering Eq. (11) for 
finite b. These artifacts would disappear for b —> oo. In Ref. 
29 they were smoothed out by averaging the result over a 
few values of b. Here, a different regularization procedure 
was adopted. The momentum sum in (11) was evaluated 
as an integral, and |%^g d | was assumed as slowly- varying 
compared to pj. Then, the integral was computed numeri- 
cally over a momentum grid and pj was integrated analyt- 
ically within each bin of the grid. 



